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ABSTR/ICT 


In "this work an a’ttemp'fc has been ijade "bo find a relatively 
simple approximate analytical solution to the Stefan’s problems 
i. e. moving interface boundary problem in heat diffusion. in 
approximation has been made in the interface boundary condition 
which facilitates to give a simple analytical solution. 

This model has been applied to first and seaend kind 
of boundary conditions • The model is applicable for a semi- inf inf 
slab. The onset of vaporisation time and the new phase zone 
width at this time has been found from the above model for 
different materials in chse of second kind of boundary condition® 
A computer program has been developed to find the numerical sol“ 
utions for the same problem without any approximation. The 
model's results have been compared with the cou^jutataonal 
results thjjs obtained and with other published results. 

Comparison shows that models's results agree with'n5% with the 
actual results for the case of meltzone width at the time of 
onset of vaporisation. But time for the onset of vaporisation 
as given by model is ns much as 50 % off from tiu o^ctual r^sT''’.ts 
for chromium. The effect of time varying heat flux on the time " 
of onset of vaporisation and oiivtti-© width of new phasezone at 
this time has been discussed. 

The computer program as developed in the case of 
semi- infinite slab has been modified for the case of cylondrica^ 
geometry. Temperature distribution, pew phasezone width and 



the time of onset of vaporisation has bean calculated for 
different materials in the case of a hollow infinite cylonder 
(infinite in length and outer radius, while with fan it e inner 
radius) with second kind of boundary condition. 



I, CHAPTER 


INTROnJClTON: 

Many practical heat transfer processes are associated 
with change of phase of material due to either melting or 
freezing, Applications include the welding of materials, 
solidification of castings and iceformation. Melting of 
solid by the application of heat flux is very iirportant phenomenon 
in case of welding. In this application the depth of melted 
zone prior to v^orisation and time for onset of vaporisation 
are the main quantities of interest. In this study the emphasis 
has been placed on these two quantities. 

The main feature of this type of problems is the 
existence of an interface which separates two regions of 
different thermophysical properties. This interface moves as 
some functicn of time. At the interface energy is either 
released or absorbed, Stefan ( 1.. ■ - ) first published his 
pioneering work on these type of problems and hence these 
problems are referred as "Stefan problems". One important 
feature of Stefan problon is the nonlinear nature of the 
interface boundary condition. Due to this# the solutions of 
diffusion equation cannot be superimposed and therefore 
particular solutions of the differential equation have to be 
obtained separately for each situation. For instance even 
if a solution satisfying the nonlinear interface conditlxjn is 
found, the same solution may not satisfy the flux boundary 
condition. To solve this problem many investigators have tried 
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approximations, which eliminate the use of diffusion equation 
in one particular phase, Goodman ( 1958) assumed the solid 
phase to be at the melting temperature and thus diffusion 
equation is to be spplied only in liquid phase while Boley ( 19 61) 
considered the case of ablation, eliminating liquid phase. But 
in both these cases all the thermophysical pix)perties of one 
phase have been igiored. In this study also the use of 
diffusion in solid phase has been avoided but the density and 
heat capacity of the solid phase has been accounted for. 

The results obtained under these assumptions have 
been compared with the anal ogue-cortqpu ter results of Cohen ( 1967) , 
Numerical solutions of the exact problem have also been obtained 
and compared with the results of Cohen ( 1967) and that of 
Hsu, Mehrabin and Chakraborty (1978). Furthermore, analytical 
results obtained from our model have been conpared with our 
numerical results* Since exact analytical solutions exist 
for the constant tenperature case, for conparison purposes, 
we have also ccnpared our analytical results with them. In 
addition, the effect of time varying flux has been discussed* 
and numerical solutions have been obtained for the case of 
cylindrical geometry with heat flux boundary condition. 

In Chapter II, literature pertinent to Stefan problem 
has been reviewed* Formulation of the exact problem and 
approximation proposed has been discussed in Chapter III, 

Method of the solution for these problems have been dismissed 
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in Chapter IV, Numerical formulation of the problem and 
problems encountered in it has also been discussed in this 
chapter. Results are discussed in Chapter V, 



II chapter 


Litei^ature Review 

The first publ iaMd' ' ow moving boundaiy 
problems in heat diffusion was that of Stefan* The original 
reference has been quoted by Carslaw and Jaeger (1). He 
obtained analytical solutions for temperature profile and the Ir.t 
interface velocity for first and second kind of boundary 
conditions* In second kind of boundary condition^ interface 
velocity was assumed constaBt. The toltial temperature was 
assumed to be phase change temperature thus making it only dne .. 
phase problem. First important exact solution was determined by 
Franz Neumann* The solutions has been given in reference C*!)* 
His solutions applied to a semi— inf in it e slab with the constant 
temperature boundary condition. Initial temperature of the 
niah was taken to be more than phase change temperature for the 
case of solidification. The solutions are in terms of error 
function and are difficult to use as a trial and error procedure 
is to be applied. 

Methods of solutions published till now are basically 
of threediff erent nature i. e* analutical, alalog computation 
and digital computation* ilv5|p G, w, (2) found the »eltzone 
width in the form of Taylor series in time. The temperature 

distribution is also given in the form of Taylor series in 

/ 

terms of time and space ilipriable* These solutions are valid for 
v%ry small time period. laiidau (3) formulated the governing 
equation in a very general form and ^rve analytical solutions fc 


effects* Goodman (4) has obtained some approximate solution 
using integral heat balance method and assuming a simpler form of 
temiberature distribution in thechanged phase. The solutions 
obtained in this paper are applical>be only when diffusion 
equation is to be applied in only one phase e*g* if material 
is initially at phase change temperature. In his future papers he 
considered some general cases also , Bol^(5) bas also solved 

the problem analytically but considering the case of ablation 
and thus diffusion equation is to be applied only in one phase. 

The most completg analysis of numerical method was 
presented by lyres in a pioneer paper* He developed lumped 
parameter system and solved on analog computer for seme~ 
infinite case with body initially at phase change temperature. 

The reference of this Paper and of Otis has been given by 
Muehlbauer, J.C.(6). Otis adopted the concept of moving heat s' " ■ 
source to account for the latent heat, Thismethod required a 
coordinate transformation in terms of an artificial time 
variable, which limits the analysis to materials initially at 
fusion temperature. After this an important work was that 
of Murray (7) in which he formulated the problem in finite 
difference form in two ways: one was of vaiaable space 
network, in which grid size varies with time to tract the 
interface position. Other uses fixed space network system, 

which has been used in this study also for numerical 
computation. Here grid size is kept constant but the boundai^- 
is tracked and the temperature of grid containing interface 
boundaiy is linearly intrapolated. The paper ^ves results 
only for a freezing problem in which the oody has been assume c 
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f 

Another important work in this field is of CohenM. I. (8) whose 
results will be discussed later in this dtudy. He solved the 
problem of melting of semi- inf inite slab due to prescribed 
constant heat flux on alalog computer. Results have been 
given for diffemt materials.. He has given results assuming 
the preperties of mat trials to be same in both phases. 

W.L. Heitz (9) has extended the variable space network 
method given by Murray so that the initialisation could be 
improved. Till new for initialisation some arbitrary new 
phase zone width was already assumed by authors. Heitz 
solved a freezing problem gnd at start he assumes all the sen=d'. 
sensible heat of subcooling is sonverted instantly to latest 
heat resulting in a solid thickness. This thicknews is then 
used to proceed further to solve diffusion equation numerically. 
C.BenacineL/ (lO) has developed three time-level implicit 
scheme, which is unconditionally stalbe and convergent on the 
basis of sn analytical approach consisting of the approximation 
of the latent heat effect by a large heat capacity over a 
small temperature range. Although apparent heat capacity 
formulations have the advantage of being simple to program, 
the predicted phase change interface location advances in 
an unphysical oscillatory fashion. Hsu, Mehrabin and 
Chakraborty (11) have solved the problem of meling due to 

laser irradication numbercally* They used variable space 
network method. They discusseii the effect of laser beam 
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intensity on the time of onset of vaporisation. Their results 
will be compared with results of this study, A brief summary 
of literature review is given in tabular form on the next 


page 



Investigator I Type of boundary { Method of Solution { Remark 

! conditions I ' 
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Otis ( 1956 ) Both Numerical ig required* Moving 

heat source assumption to 
account for latent heat 
effect* 





III. .CHAPTER 


Exact formulation for the Stefan problem of semi- 
infinite slab with different boundary conditions and for 
infinite hollow cylinder with heat flux boundary condition 
is discussed below, 

3,1 Semi-infinite Slab 

Consider a solid semi— infinite slab initially at a 
uniform temperature which is assumed to be less than the 

melting ten^erature The geometry and the coordinate system 

used ais shown in Figure 1 ( a) R heat flux Q( t) is applied to 
its surface at x = 0, 

The governing heat diffusion equation and boiandary 
and initial conditions as given in Ref, ( 4 ) can be writtai as 


r>T(x,t) 

- cy ^ 

“ 2 

f X 50 

(1) 

~?T( 0 ,t) 

.>X 

= Q(t) 


( 2) 

V co.t) 

>)X 

= 0 

ft 

V 

0 

C 3 ) 

TCx, 0 ) = 

^o 

; X >0 

( 4 ) 

"3 T( X , t) 

;y 

J x( t) >X70, ty7t>tj^ 

( 5 ) 

^ ~gTCx.t) 
^ 7 ?X 

== Q(t) 


(6) 
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T (X( t),t) = 




>t 


>t_ 


m 


~a-T(x^t) 

■at 



>2 


T(x,t) 




X>x(t), V>t>tj^ 


with interface condition 

iLJSCxTc t),t) ;^x+( t),t) dX( t) 

”* sa?x dt 



(7) 

( 8 ) 


C9) 


The system of equations (l)-(9) fully describe the 
problem of melting of semi— infinite slab due to prescribed 
heat flux* In case of constant temperature boundary condition 
eguations ( 2) and ( 6) will change to 


T(o, t) = , t > 0 

S 


( 10 ) 


3,2 Hollow Infinite Cylinder 


Consider a hollow infinite solid cylinder of internal 
radius Y'l. at a uniform tei^erature which is assumed to be 
less than melting temperature The geometry and coordinate 

system used has been shown in Fig,l(d), There is only the 
difference of gecmetry between the previous problem and in 
this one. The governing equations and boundary and initial 
conditions can be written in a similar way as follows: 


3*' T( r^t) 
^t 




1 2>T(r.t) , 

r > r 




( 11 ) 



3^ T ('i2#t) 


T( OD , t) 

s ^ r 


Q(t) 7, t^'>t>0 


( 12 ) 
( 13) 


= 0 


7 t >0 
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T( r,0) S T 


o 




( 14) 


1 ?>T( r.t) 




m 


- ^ 


'^Jli r . , t) 
^r 




( 15) 
( 16) 


T(R(t),t) = T 


m 




m 


( 17) 


^ r.t) ^ ^ i ( 18 , 

O' IT 1C 


with interface boimdary condition 


K 


*(t),t) _ >T(R+(t),t) . o, dR(t) 

- “'^s ^ +/ ^m “"dt" 


( 19) 


^inalytical solution of equations (1)— (9) is difficult 
to obtain because cf the interdependence of diffusion equations 
(5) and (8) through a nonlinear boundary condition (9), Problem 
becomes more complicated due to the existence of heat flux 
boundary condition at the surface, Evans G.w, has given 
analytical solution to this problem in the form of Taylor 
series, but solution can not be applied for large times as the 
series diverges. Other att^ipts e.g, by Goodman and Boley 
simplify the problem to one phase 1^ suitable assumptions 
which are discussed in Chapter II, In all these attempts all 
the properties of the solid region have been ignored . in 
this thesis an attempt has been made bo decouple the diffusion 
equation in solid region from the rest of the system. The 
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interface condition (9) is modified such that density and 
specific heat of solid is taken into account tut the diffusion 
equation is not required to obtain melt zone width and the 
time for onset of vaporisation as discussed below. 


3,3 Proposed Model 

Consider the interface boundary ccnditicn ( 9) , The 
left hand term denotes the heat flux caning fron the liquid 
region and the right hand side terms denote heating and melting 
of solid. The proposed model is obtained as follows. Consider 
the grid system as shovjn in Fig, 1(b), Just at the. onset of 
melting we can write 


-K 


7 > T 
S ,'^X 

■o 




dX I 
dt 


K 


S '^X 


C 20) 


= c (T-TO 


dX 


r-p ' -"m "I'dt 




-K 1 

s ^x ^ 


2^ 

C 21) 


in a similar way 




^ / m i-r* C3,j Xr f 






V ' "■n+2>i| 




- K . 

S -=^x 1 

{ 22) 



CD 


but 

OT 

s= 

s ^ 

0 and T(cq t) = T^ 

( 23) 
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Hence assuming is sane for all grid points we can write 


(X'^(t),t) = c (T - T ) 


s 


p m o 


Thus the condition (19) changes to 

(3r(t>,t) =po(T .T)^^PL 

^ J p m o dt J m dt 


K. iL2!CX“(t),t) _Q * dX(t) 
^ '^x "J dF 


( 24) 


( 25) 


(26) 



L^, henceforth will be referred as modified latent 
heat of fusion. 



IV, CHM>TER 


Method of Solution 

In this study the Stefan problOT is solved in order 

to find the melt zone width and time for onset of vaporisation 

for different materials. Therefore calculations have been 
2 

made to find , and Al X,,, where AI is the absorbed flux 

and'-’J^V ^ time between onset of melting and vaporisation 
and melt zone width at the onset of vaporisation respectively. 
These values have been foiand to be constant for a particular 
material. In the calculations made here, the properties of 
the materials have been averaged out in the temperature range 
ccncemed e,g. in case of solid the temperature concerned is 
from T^ to T^ while for liquid it is T^ to T^; 


Model 


To find the analytical solution of the problem of 
semi- infinite slab with heat flux boundary condition equation 
(l)-(8) are to be solved with the modified boundary condition 
( 26) , This problem is now broken in two segmCTits. One upbo 
t = tjj^, of which the standard solution as given in Ref, (4) is 


T(x,t) = I - y erfc y) 

s ^ rfH 

where 


y = x/( ? erfc y = ■ 




/'CO ^ 

] e-2 dz 

Jy 


( 28) 


( 29 ) 


Here Q is assumed to be a ccnstant heat flux. 
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After the onset of melting two phase region exists. 

To simplify the solution without much loss of accuracy a 
parabolic distribution is assumed in the liquid region as 
follows: 

T(x^t) - A(x-X{t))+ B( X - X( t) 4 - Tm (30) 

Thus now the solution obtained by Goodman can be 
applied here. The results are being quoted here from Ref ,( 1) , 


where 


and 


A = ifl - (1 + 4 


1/2 


I/O 1 


L 

m 

X( t) 
pL 


B = i- 


8 


Ju.= 

I 


Q(t) X(t) 

p^m 


( 31) 

( 32) 

( 33) 


Q(t) lQ(t)dt ^r 

“*■ = 'e fM+ 5 + 


) m 


( 1 + (34) 


Surface temperature can be obtained from ( 30) by 
putting X = O and values of A and B, This gives 


T 





(35) 


Thus the solution obtained by Goodman for the problem 
of initial temperature T^ = T^^ can be applied even if T^ is 
less than T^^, by simply modifying the latent heat of diffusion 
as shown in eq.(27). 
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To verify the validity of the modified latent 
heat assumption, this simplification has also been applied 
for the case of constant surface temperature. Exact analytical 
solution as given by Neumann are being quoted here for coitp arisen 
from Ref, (4). 


C T - T ) 

T(x,t) = Tg - — erf ( 


X 


erf^c 




-)? X(t)>x>0 


T(x,t) = 


erfc 


erfc ( — ); x^X(t) 


c'sfer. 


where /\, is given by 


-A? 


(Ts - T^)e-’^ C Kg 

erf ^ S 


\ 2fC: 

)e'’^ X 


er£c(A 


2f^ 

® -h HnAc 

+A "c 

pL 


( 36) 

(37) 


( 38) 


and interface position is given by 


X(t) = 2A^ eXVt)' 


(39) 


While v;ith the use of the proposed model, equation 
(9) changes to equation (26) and thus the solution can be 

X 

obtained in a similar way as obtained above by Neumann, In this 
case equation ( 38) will change to 


(T^ - T„)e 
' s m 


-A 


Ll X 


Af 


= n 


m ' •■c 


' C 


(40) 


'PL 


which is in much simpler form, other equations will remain 
same except the value of changes to obtained from (40) 
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Numerical Scheme 

To solve the exact problem given by the set of 
equations (l)-(9) numerical oDmputation has been carried out, 
so that the analytical results obtained frm the model can be 
compared* Three point central difference formulation has be^ 
used in this scheme so that algebraic equations obtained^ 
due to the implicit time schone used, are in tridiagonal form 
and are much easier to solve. Figure 1(b) and ( c) shows tte 
fixed space network used. This procedure has been taken from 
Murray but he used parabolic interpolation to track the interface 
while here for the sack of simplicity without much loss of 
accuracy linear interpolation has been done, other method 
given by Murray is of variable space network. In this tracking 
of interface is not difficult but the size of the grid varies 
with time. Implicit technique has been used to avoid the 
limitation on time or space grid size. In case of explicit 
method which is easy to ccnpute this limitation can not be 
avoided because otherwise instability will arise in the 
solution. 

First of all the variables have been non-dimensicualised 


as follows: 

Subscripts being dropped where they are obvious from 
the context. 


O 



■k 


X 


xA ; t 



(41) 


diffusion equation becomes 

t* =3 Ax* ^ 


( 42 ) 
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in finite difference fomi 


"^^i-l + + 2 ) 1 ) 


= «i 


(43) 


where 


•\ ^ 

A = At /ax 


*2 


(44) 


The space grids have been chosen as shown in Fig.lCb), 
Since first grid has been taken of full grid size the botindary 
condition (2) changes to 

Q(t) = ^ |S+ C)qpAx ^ (45) 


in non-dimensional form it can be written as 


c?' t 


Q( t)L 

K(T.^-T„) 
V m 


and in finite difference form 


(1 +/\ ) ©T — A 


^©2 = «1 + 


Q( t)L ^ t 
* 

Ax K( T^-T^) 


other boundary conditions ( 3) and ( 7) can be written as 


(46) 


(47) 


®N-1 “ ° 

© (X(t)) = o 


(48) 

(49) 


The moving interface boundary conditicn (9) in non— dimensional 


form becomes 

%■ 

0X 


X=iX"” ( t) 


■'s 


(50) 


otdC'^(t) 


where 
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V = vAe ? = 


and 


V s 


(3X( t) 
dt 


( T - 


m 






( 51) 


(52) 


In this scheme fixed space network, method has been 
used hence it is not necessary . that interface always lies on 
same grid point. Thus the interface is to be tracked at each 
time step and the temperature of the grid ccntaining it is 
interpolated linearly. Thus if x is the distance of the 
interface from the grid point then equation (43) can be written 
in finite difference form as follows 


it 

V (t) 


K 

K, 


...■0(nr + 1) 

* T 

(Ax - Sx ) 


Q inr -» 17 

4“ ★ ★ 

(Ax ) 


( 53) 


Ax 


^^6x 


In cylindrical geometry case the definition of 
variables changes as 

^ t 


* V * 

r = r/( - r A ; t = 


( 54) 


(r„ - 


Diffusion equation (11) in nondimensional form is given as 


'd 


Q 


2^0 
?2 


5 Q 


(55) 


dt dr " r Or 

which in finite difference form reduces to 

“^i-1 ^ ^ + 2)\+/\ ^ 


= O. 


( 56) 


A =AtVAr*^ 


where 


(57) 
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Rest of the formation is same as in case of semi, 
infinite slab. 


Initialisation and Tracking of Interface 


In the case of numerical computation till the onset 
of melting there is no problen in computing temperature 
distribution for differ^t time steps. As soon as the surface 
reaches the melting temperature, two phase problem is to be ' 
considered. Since at the start there is no liquid zone, the 
diffusion equation can not be applied there. Therefore first 
the velocity of interface, which is initially at x = O is found 
from the boundary condition (9), it is given as 


because 


v(t) = 


Q C t) 

7^ * vi 


■3k t) 

X 


X=X'^C t) 


T 

~ at X = O 


( 58) 


( 59) 


This equation is applied till the first grid is 
covered by liquid region. After this the temperature of the 
first grid point is found by taking energy balance for it, as 
given by equation (47), The tempera "ture of the second grid 
point is interpolated between the first grid point and the 
interface temperature if it is in liquid region and between 


the third grid point and interface temperature if it is in solid 
region. This procedure is repeated till three grid points 
are covered by liquid region. 


equation can be applied in liqui 



Aoc 


.No. 


ds three points, 
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The temperature of last point in the liquid regicn and first 
point in the solid region is found hy linear interpolation 
between interface and the previous or next point respectively. 
Main difficulty in the numerical computation comes 
in calculating the velocity of interface from eq,(53). 

This is because if the time step is large and the grid size 
is very small then in a single iteration the interface may 
cross more than one grid point. Since the temperature of the 
solid region was quite low and if we calculate slop in liquid 
region at boundary fron this position it will come quite low 
and thus may result in negative velocity. Consequently, the 
time step is to be reduced or the grid size is to be increased. 
Thus even thou^ the numerical scheme used here is in implicit 
form there is a limit on time or grid size to avoid this 
unphysical oscillatry motion of interface boundary. Moreover 
chosing. Ax larger may result in higher error in computation. 


Error Analysis 

There are two types of errors in numerical computation 
C i) Error due to finite difference formulation 

( ii) Truncation error 

In finite difference formulation of the diffusion 

2 

equation the error will be of the order of o( /At) + 0( Ax) 
which for the. present calculations comes out to be less than 2%. 
Truncation error arise due to the finite number of significant, 
digit calculations. Since in the calculation ei^t significant 
digits have been taken by the computer that itself is quite 
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accurate and if there is any error at the last digit place 
it is almost insignificant in comparison to the error made 
due to finite difference formulation* 



V, CHJ^TER 


RESULTS AND DISCUSSION 

First we discuss the results for the case of semi- 
infinite slab. Figs, (2) and (3) and Table ( 1-a) show a 
compariscn between our analytical results, numerical results 
and the results of Cohen, Since Cohen obtained these results 
assuming same properties in both phases, for comparison 
purposes in Figs. ( 2) and (3) and Table ( l-a) we have also 
done the same. It can be seen from these figures and table 
that the values of AI^X^ and AI^'(v are constant for a given 
material. This result is well known and was also obtained by 
Cohen, From Fig,( 2) we see that values of AI X„ obtained from 
our analytical model agree within 5% with those of Cohen and 
our numerical results, Ne also see that the error in model 
results is positive for some materials and negative for others. 
Thus there is no definite trend and it is felt that the 
difference in results is so small that the positive and negative 
errors may occur due to inherent error in finite difference 
formulation. 

The values of obtained from our model differ 

as much as by 30% ( for instance in the case of chrcanium) 

from the numerical values. For most of the materials 

discussed this difference is about 20%, Error percentage in 

2 ■ 

the values of snd Al^^y are given in Table ( 1-b) for 

different materials. From this table we see that the error 
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percental in the case of is generally higher for those 

materials which have smaller values of thermal diffusivity, but 
some exceptions to this effect can also be seen. For instance, 
in the case of silver though diffusivity value is comparatively 
large, the error is also large. Similarly in the case of 
Platinum diffusivity value is comparatively low but the error 
is also low. 

We have also obtained results for those situations in 

which the thermal properties of two phases are different. These 

results are given in Figs.(4) and (5) and Table (2), It can 

be seen that these results differ significantly (as much as 

by 155% for Gold) from those obtained by assuming constant 

properties in both phases. Therefore to obtain realistic results 

properties 

it is necessary to take into account different ,dn different 
phases. 

Hsu et»al, have also obtained numerical results 
considering differoit properties in the two phases. In Fig, (4) 
we compare our numerical results with these of Hsu et.al. 

Even though Hsu et,al, used variable space network scheme while 
we used constant space network sdhone, within 3% the two results 
are the same* Figs, ( 6) and (7) give the variation of and'^^y 

with respect to absorbed flux for different materials. The 

—1 -2 

results show that x -<*->1 and T ’- rr ^ I ^ as discussed above, 

y O V o 

In Figs,(8) and (9) we show the effect of the flux-pulse- 
shape variatim on '^y and Xy, Two shapes one r^tangular 
and the other Gaussian have been considered. Here the maximum 
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intensity for the two cases has been taken as In 

comparison to the flat topped pulse, the Gaussian pulse takes 

more time to get the material upto the vaporisation stage. 

Further the molt zone width is also more with the Gaussian 

■pulse. But looking from the en'ergy point of view the rectangular 

2 

pulse takes 285,76 J /can while the Gaussian pulse takes 
2 

744,36 J/an to reach upto the vaporisation stage. This is 
2,6 times more than in the case of flat topped pulse, while 
the width melted is^ only 2,22 times more. Thus we see that if 
the quantity of interest is the melt zone width per unit of 
available energy than the rectangular pulse is better. 

Table (3) gives the values of \ for the case of 

' "C 

constant temperature boundary condition. Both these results 
have been obtained analytically, one with the proposed 
approximation and the other without any approximation as 
suggested by Neumann, The difference is more than 50% in most 
of the cases which shows that the model is not applicable for 
this case, ■ One reason for this difference is that in this 
case there is no concept of molted thickness at the time of 
onset of vaporisation. 

We have also obtained results for the case of cylindrical 
geometry numerically. These results are ^own in Pig.{ 10) , 

From this figure it can be seen that even for r^^ as small as 
0,15 cm the values are very near to that of semi- inf inite 
slab case for Aliaminium, This shows that analytical results 
obtained for the case of semi- inf inite slab can be applied here 
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if the internal radiuns ig. relatively large, in general the 
numerical calculations show that if the internal radius is 
more that approximately ten times the melt zone width then 
the values obtained are very close to that of semi— inf inite 
slab case, 

CONCLUSIONS 

In this study an attanpt has been made to find a 
relatively simple analytical solution to Stefan problems i.e, 

i 

moving interface problems in heat diffusion. An approximation 

was made in interface boundary condition to delink the liquid 

region from the solid region aich that heat capacity and density 

of the solid phase is still accounted for. In the case of 

constant heat flux boundary condition the melt zone width 

at the time of onset of vaporisation obtained from the model 

is found to be in good agreement with the exact results 

obtained numerically and analog ctxnputer results obtained 

by Cohen, Percentage error in the values of Al^^Xy is less 

2 

than 5% but in the case of it is sometimes as high 

as 30% as in the case of Chromium, But generally this error is 
around 20% for most of the other metals. In the case of 
constant tenperature boundary condition, however, results 
obtained frcm the model differ from the actual analytical 
results by more than 50%, Numerical results have been obtained 
for the qase of qylindrical geometry with constant heat flux 
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boundary condition. Results show that if the inner radius of 
the hollow cylinder is more than approximately ten times the 
melt zone width then the results of semi— inf inite slab case 
can be applied. Concerning the improvement to the above model 
it is desirable that the thermal conductivity of the solid 
is taken into account as this property has not been included in 
the above model. 
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semi- in finite solid with constant heat flux boundary condition 
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TABLE 3 


For the case of semi- inf Ini te solid with cx^nstant 
heat flux boundary condition 


it 

Materials 

X 10-2 

T / % 2 

J /sec on 

^ ■ 

1 

\% difference 

1 -W ^ - 

^ i 

iResults of Ref, No, (11)1 

Nr ,C ompu ta tion al 1 two results 
results .* 

Aluminium 

4,21 

4.33 

2.85 

Iron 

4,07 

4,02 

1.24 

Nickel 

4.27 

4.16 

2.644 

* 

P roper ties 

of materials taken are 

given in Table 

( 6), 
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imiE 4 

For the case of semi- infinite solid with constant 
surface temperature (=Tg) boundary condition 


* 

Materials 

1 

Model* s result 

X** 

1 

1 

t 

{ Exact analytical results 

i 7\i* 

Chromium 

0,895 

0,53 

Copper 

0,910 

0,63 

Gold 

1.01 

0.74 

Iron 

0*952 

0,62 

Nickel 

0,77 

0,50 

Platinum 

1,01 

0,70 

Rhenium 

1,045 

< 0,70 

Silver 

0,96 

0,655 

Tantalum 

0,837 

0,53 

Tun gs ton 

0,888 

0,56 

Aluminium 

1,09 

0,81 


Properties of materials are taken as giv^ in Table (5)« 

ii||^ ^ 

^ is non*- dimensional variable given as 

X( t) 
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t=Cu Assuming block body and 
2=:Alf some properties jn both 
3= CrJ phases. 

4 = Al with diff. properties in two 
phases and 
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Tnls prograiii computes trja temperature aistrioution, tne 
posltioni of Interface and interface velocity wltn resoect 
to tlme.iThe body is one-dlmenstonal finite or settl „tntf nite 
Slab at a constant temoeratare TO Initially. arescrfael 
beat flax Is applied to its one face and otner bsino insulated, 
rns program is divided in three portions, First ons calculates 
fne temps, dlstrioution in the slab u.ntill tne on, set of 
msltinig. Second one taices it as initial temo, dlstrioution and 
computes tamp, distribution, interface position and Its velocity. 
To solve the heat conduction eouation finite difference form 
with central diff,' has been used, Third portion is a subroutine 
to solve; simultaneous linear algebraic equations resulting 
from finite difference formulation <flth central difference. 

The material properties except density has been ta<-en diff, 
in twp phases but in a sinale ohase they are const, '•^Itn 
temp,' 


hIST DF 

MasAppiiedl heat flux.' rgslnltlal temoerature.' 

T'lsHeltihg temp, tSa.VaPorlsation temo. 

C^HStCOiih are thermal conductivities in solid £ lid,- reso, 
CPS.CPli are heat capacities in t«o phases, 

MS.ABSif' are flux aOsorPtlvltles in two phases, 

DIFS, Dirt* are thermal dif fusl vltles In two ohases. 

DfiHS is density, tatcen same for both ohases, 

It'fshateht hear of fusloh,' 

M*«Slab width, 

AlMAXsMaximum flux in case of time varying flux, 

BDsdalf width at half max.yflux variation is Sausslan, 

Program has been written In non-dimensional form,' 

Properties values are in ".3,5, unit system, Temos, are in degree 
Kelvin.* 
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08600 

08700 

08800 

08900 

09000 

09100 

09200 

09300 

09400 

09500 

09600 

09700 

09800 

09900 

10000 

10100 

10200 

10300 

10400 


10 


20 


30 


DlMEOSIDfi TA(IOUO) , T ( 1 000 J , A C 1 000 ) , B( 1 000 ) , C( 1 000 ) , D C 1 000 3 , 
9KXtl0003 ,DA{10003 

A10 = 1.0K+05;AIMAX=1 ,OE:+07|BO=0.0005?TIMD=BOtSORT(2.03 
TU = 300.0,*AB6=1,u;ABSL=1.0 

READ*,DENS,COwS,CPS,CONb,CPL,XM,TB,ALM,XI 

N = 1000?iM=1000 

!'(M = i'>i-l ; NMjM = in-2? I^P=Nfl 


initialisation and nondlmensionaiisat ion of variables. 


TlMt;;=0.0?DlST=0.0;DT=0. 000001 

IfeRl=0;!UT=0?IiMTl=30;lNT2=10 

ITER2=0 

A6=1.0 

A8S=5,03 

CONL.=COt>iS? ABSb=ABS;DTFL=:DIFS?CPb=CPS 

A8SB=0.15 

CUNSL=CONS/COlML 

DX=Ali/N 

DXH=0X=<'3. 0/2.0 
DXT=OX*5. 0/2,0 
DXfl=DX/AL 

OIFS=CONS/CDENS*CPS) 

OXFL=COML/CL>EiMS*CPL) 

DrMS=DirS*OT/(AL*AL) 

D rJ'IL=l>IFL»DT/ ( AL»AL) 

A6AMS=DrNS/(DXN^DXN) 

A6AML=DTNL/CDXN=I=DXn5 

ALMMsALM+CPS*(T.'>Tn) 

VCM=C0N6=t'(rb-Tf1)/(DENS»ALMM*AL) 

VC=CONL+(TB-T^i)/(UENS»ALM*AL) 

Y=ALM/(CPS»(TM-TO)3 

PRIMT10,M,AI0,AL,DX,DT,Y 

F08MAT(5X,' N= ' ,I5.5Xf 'AIDS SE15,8,5X#' AL® 
8= %F7.4,5X,' DT= *,£15.8,' Y= ',F5,2,//) 


,F4.1,5X,' DX 


Solving ciitfusion eq, till tne onset of leiting. 

DJ20 3=1, j 
ACT 3=0.0 
HCI3=0.0 
C( I )=U,U 
0113=0,0 
CJ'.n'lr.UE 
B(U = 1 .OtALAfIS 
C(l) = -ADA.-iS 
0030 1=2, 

Aa3=-AL.A^5 
8CI)=:1.0+2.U + A3.AiVS 
C(I) = -Al4Ai-lS 
At'13 = 1.0 
tJtNjs-l.O 
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11200 
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W 

11300 
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11400 

Iv* 

11500 


11600 
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11700 

45 

11800 


11900 


12000 

50 

12100 


12200 


12300 


12400 


12500 

90 

12600 

100 

12700 


12800 


12900 


13000 


13100 


13200 


13300 

101 

13400 


13500 

110 

13600 


13700 

120 

13800 


13900 

130 

14000 


14100 

132 

14200 


14250 

r» 

w 

14275 

L, 

142H7 

w 

14293 

/% 

143U0 


14400 


14500 


14600 


14700 


14800 


14900 

•H 

V, 

15000 

W 

15100 

r% 

15200 


15300 

135 

15400 


15500 


15600 



D04O 1 = 1, N 

r(I.) = (TO-rM)/CTB-lM) 

Din 10 jj=i;,M 

TlHE=IIMEtDT 

JJMsJJ-1 


If flux Is bell stiaped.AIU is given as, 
A£0 =AIHAX=»EaP("C lTIME-TIMO)/8U5=f*2) 


C0=UTi'jS'f AbS!fAlU<'Ali/( CTB-T«l)*COflS*UXN) 

PRIiyT45,CU 

F0RMAr(5X,E15.8,/) 

OCU=CO+TU) 

00501=2, il«l 
0(1)=T(I) 

IfERl=irEHl +1 

CALL THXDAGC A,B,C,D,N,XX) 

D09O 1=1, W 
TCI) =XXC I ) 

T A C 1 ) =T £ I ) ♦ ( TB-TM ) +TM 
FORMAT (5X, F8. 1 ,10X,fc;i5.8) 

MIT=NIT+r 

TS=(3,0»TCl)-.TC2))/2.0 
IFCMIT.LT.IMTll GO t6 101 


MIT=0 

TSA=TS»CTB-TM)+TM 

PRlNriOO,XSA,TIME 

CONTINUE 

irCXS.GE.O.O) GO TO 130 

CONTINUE 

PRINT120 

FORMATClOX, 'Temp not reached to melting 
STOP 

CONTINUE 

PRINT132 

F0RMATC5X,' TIME TS 

921, DIST *,//) 


//) 


TN 


Aelting nas started. Calculations for liquid zone growth till 
Tnree grid points are covered. 


i)r=uT*2.0 

DlNSsOirSfOT/CAL^^ALJ 

0T,Mj=DIFLfi)r/{4ij=fAL) 

A L A’A 5=0 r N S / { n X N * OX N ) 

ALA*iL=U'i. Uj/(DX ' i^fDXN) 
CS=DT=fA0Si‘A10/(DEiiSfCPS*DX) 

CS.« = UT:'i3=f AHSfA£0+AL/((Ta-TM)*CONS»DXN) 
Cu=Ur^Ab5Li=A I.U/(DENS*CPL*UX) 

C-j-JsDrNL + AaSL^Ain^AL/t ( 10-T'O^CONL»DXN) 
Pnl 't n 35 . T 1 MF 

FJR'-'.ATCSX, ' rime for melting ',E15.8,/) 
AiO=AIMAX*EXPC-(CTIi4E-TiMO)/B03**2) 

C Jln = A.35inAin*AL/(CONjU*(TB-i’M)) 
v/=Cai.n2,OfCON5u^(Tll)-TS)/DXN 


lb7U0 
15800 
15900 
16000 
16100 
162U0 
16300 
16400 
165U0 
16600 
16700 
16800 
16900 
17000 
17100 
17200 
17 300 
17400 
17500 
17600 
17700 
17800 
17900 
18000 
18100 
18200 
18300 
18400 
18500 
18600 
18700 
18800 
18900 
19000 
19100 
19200 
19300 
19400 
19500 
19600 
19700 
198U0 
19900 
20000 
20100 
20200 
20300 
20400 
20500 
20600 
207O0 
20800 
20900 
21000 


140 


150 


160 



rf 

i0 


c 

r* 


IbS 

1 70 
180 


7A=V*VC 

VA=ABS^=AUJ/CD!':^'IS*ALM^U 

D1ST=D1ST+7A*UT 

OlST?4=DlSr/A0 

IF(DIST.GT,DX/2.0) GO TO 170 

DeijX=:DX/2,U-D{ST 

oaox?i=Dfe:LX/Ai, 

A(i)=0,0 

8(1) = 1.U+1jX/DE1>X 


CCU=-1.0 
5 - 0.0 


0(1. 

00150 1=2, UH 
A(I)=-ALAMS 
rt(l5 = l,0 + 2,01'ALAMS 
C(I)="ALAi4S 
0(U='i(l) 

ACMlsl.O 
8CM)=-1 ,0 

ct’O= 0.0 

0(N)=O.U 

TIME=TIME+DT 


irEH2=riEK2 + l;faT=NXT+l 
CALL TRIDAGC A,B,C,D,N,XX) 
00160 1 = 1, 

TfI)=XXflJ 

TM1)=TCI}*CT8-TM)+TM 
IF(NIT.LT,IiyT2) GO TO l6l 
NIT=0 


PRINT100,TA(1),TIME 

COMTIOUE 


AID as oelow in case of bell shaped enercty pulse, 

A10=AlMAX^EXPC-C(TlME-TIMU)/80)**2) 

COIN=A8SL*AIO*AL/(COf4L*tTB-T?4)) 

V=COIN+CONSL*TC 1 ) /DELXN 
VA=V*VC 


Interface velocity In model's case is as oelow, 
VA=ABS*AIO/(U£NS»Af.MM) 


0IST=LIST+VA*1>T 
ISA=(3.01‘TA(1)-TA(2) )/2.0 
PRIrjT165,TlWK,T6A,TA(i'l) ,VA,DIST 

FJRmAI( 3X,F:15.8,3X,F7.1 , 5X , F7 . 1 ,5X , fcl5 . 8 , 3X , E15 .8 , // ) 
IF(H3T.GI.L/X/2.0) go TO 170 
GJ TU 140 

IFCUIST.GI.DXH) GO TO 210 
i)lST?I = UiST/AL 
DEi,X=DlST-DX/2.0 
nXLXO=OELX/AL 
TS=UlSTiaC01iJ 
A CD =0.0 
BCU = 1.0 
C Cl 3=0.0 

DCUsDELXiaTS/DlSTN 


21100 


A(2)=0,U 

21200 


8 ( 2 J — 2 « 0 "DEL X/ OX 

21300 


C(2j=UEL.X/DX-1.0 

214UO 


i)(2) = 0,a 

21500 


OU190 1=3, OM 

21600 


ACI1=-ALAMS 

21700 


B(l)=l .0+2.04ALAHS 

21800 


CU)=-AOAMS 

21900 

190 

D(1)=T(1) 

22000 


flME=TIME+DT 

22100 


lXEH2=IIEH2fl;Nll=(«lTf 1 

22200 


CALu THI0AG{A,8,C,D,N,XX) 

22300 


L)U200 1 = 1, N 

22400 


fCl)=XXCIJ 

22500 

200 

TA(l)=TCIH(XH-TMl + rM 

22600 

s-> 

PRIiVnUO, CTACl) ,I = liN) 

22700 


IFCiMlT.0T.iKT2) GU fo 201 

22800 


NlT=0 

22900 

C 

PHXNriOO,TA(l) ,T1ME 

23000 

201 

CONTINUE 

23100 


V=rCU/UELXK+C 0 KSL*T( 2 )/(UXM-DE 0 XfO 

23200 


VA=V*VC 

23300 

C 

V=TC1)/UELXK 

23400 

c 

VA=V*VCm 

23500 

c 

PRIKT165, VA.TXME 

23600 


U1ST=01ST+VA1'0T 

23700 


rSA=(3.0»TA(l J-TA(2))/2,0 

23800 

c 

PRINT165,XXME,T6A,TAIN),VA,01ST 

23900 


IFCDIST.GT.0XH) GU TO 210 

24000 


GO TO 180 

24100 

210 

IFCOXST.GX.DXT) GU TO 240 

24200 


DeLX=DISr-DXH 

24300 


DELxN=t)ELX/AO 

24400 


A(l)=1.0 

24500 


BC1) = 1.0+ALAML. 

24600 


CC1)=-ALAmL 

24700 

m 

A10=AlMAX^EXP(-( ( •XIME-TIMU)/B0)=i‘^2) 

24800 


C OMsOTMLi^AdiiL^ AlU^ AO/ 1 (XP-TM) ♦CQNL=<‘UXN) 

24900 


0(1 ) = 1(1 )+CLAi 

250U0 


A(2)=1.U 

25100 


8C2)=-1 ,U-DX/OEL.X 

25200 


C(2)=y,0 

25300 


o(,!:)=u,o 

25400 


AC3)=0.0 

25500 


P(3J=2.0-uEl.X/!)X 

2560U 


Cti)=OELX/uX-l.U 

257g0 


i.>(i)=u ,u 

25800 


1)0220 1 = 4, svw 

25900 


AC n=-Al.A,-lS 

26000 


8Cl) = i.0 + 2.U*AI.APS 

26100 


CCI)=-ALA.‘(S 

262UU 

220 

OCX J=XCi) 

25300 


ri*'lh=XiMK+UJ 

26400 


I rEH2 = IXt;k2 + l;MlT=KiT+l 

26500 


call. ‘lKiL'AGCA,b,C:,D,N,XX) 



26600 


00230 1=1. N 

267O0 


r(IJs:XX(I) 

26800 

230 

TAtl)=T(l5^(TB-TM)+rM 

26900 

27000 

C 

PRINTIOO, CTA(I) ,I=1,N) 

IF(iUT,6T.lNT2) GO fo 231 

27100 


f<lT=0 

27200 

W 

PRINT100,TAC1) ,TIME 

27300 

231 

CONTINUE 

27400 


V=Tt2)/DE’lLXN+COWSL^T(3)/(DXN-DeLXN) 

27500 


VA=V*VC 

27600 

C 

V=T(2)/DE0XN 

27700 

C 

tfA=V*VC« 

27800 

27900 

28000 

c 

PRlNTi65,¥A,TlME 

01ST=0ISTtVA^DT 

XS A= C 3 . 0 A ( 1 ) -TA 1 2 ) ) /2 , 0 

PRIM T165/riME,'ri>A,TA(N),VA,D 1ST 

IFCDIST.GI.OXT) GO TO 240 

28100 

28200 

n 

C 

28300 


GO TO 210 

28400 

240 

CONTINUE 

28500 

242 

TIMM=( (TM-T0)/C2,*ABS»AI0) )*’*^2»3.i4*DENS*CPS=l'C0wS 

28600 

C 

TIMM gives analytical value of tm , 

28700 


PRINT260,TIMM 

28800 

260 

FORMATCSX,' SE15.8,//J 

28900 

C 


29000 

c 

Atleast tnree griti points nave come in liquid region 

29100 

c 

Calculations nereafter. 

29200 

29300 

c 

NJJsNOJ+l 

29400 


00315 JJ=NJJ,M 

29500 


JJMsJJ-1 

29600 

c 

PRINT280,D1ST 

29700 

280 

FQRMATCSX,' DI5T= ',E15,8,//) 

29800 


Q=OIST/DX 

29900 


D0290 1=1, N 

30000 


NR=1 

30100 


NRMM=NR-2 

30200 


NR«=NR-i 

30300 


NRP=NR+1 

30400 


NRPP=rtH+2 

3U5U0 


IFCO.LE.l) GO TU 300 

30600 

290 

CJNTii.UE 

3U7UU 


GO TO 500 

30800 

3U0 

i>fcLX=UlST-Ci^R=»‘UX-UX/2,0) 

30900 


OEBXNsD&BX/AG 

31000 

V 

P H I N T 3 i 0 , »* K 

31100 

310 

FJRi^AUSX,' NK= ',I10,/J 

31200 


00120 I=i,mR 

31300 


A Cl 1=0,0 

31400 


B(I)=0.0 

31500 


cti J=O.U 

31600 

31700 

320 

DC!) =0 . u 

C 0 ‘i T 1 w U £ 

BCU = 1.0 + ABAML 

31800 


31900 


CC U=-ALAMU 

32000 

r* 

AIO=Ai.'aX^£XPC-(C 11HE-TIM0)/B0)**2) 

32100 


C6N = l)T NL»ABSL*AlU»Al,/t t TB-TMl+CONB^DXN) 

32200 


DCU=CBN + TCU 



32300 


It''’(i)t’.6X.bT.0.0) GU ft) 350 

32400 


03330 1 = 2 ^ N 8 M 

32500 


A(IJ=- ALAMO 

32600 


B(I)=1 .0+2.0»ALAML 

32700 


CCI)=-ALAHL 

32800 

330 


32900 


ACNf<) = l ,0 

330U0 


6UR)=-i ,0 -DXN/IjELX« 

33100 


CCNK)=0.0 

33200 


D(NR)=0.0 

33300 


ACNHP)=0,0 

33400 


l3tNP.PJ=2,O-0ELX/DX 

33500 


CCNHP)=DELX/DX-1 ,0 

33600 


DC. MRP 5 = 0,0 

33700 


D0340 I = NRPP,IMH 

33800 


ACr)=-ALAMS 

33900 


BCl J=1 ,0+2,0*ALAMS 

34000 


CCI)=-ALAMS 

34100 

340 

OCIJslCi) 

34200 


A(rJ) = i.O 

34300 


BCN)=-i,y 

34400 


CCM)=0,0 

34500 


DC«)=0.0 

34600 


GO TO 380 

34700 

350 

00360 I=2,RRMH 

34800 


ACI)=-ALAML 

34900 


B(I)=1,0+2,0*ALAML 

35000 


C(I) = -ALA.'IL 

35100 

360 

oa)=TC13 

35200 


A(NRM)=1 .0 

35300 


8(NRM)=-C2.04:dX+D£LX3/(DX+DELX) 

35400 


CCNRM)=0.0 

35500 


D£N[RM)=0.0 

35600 


A^MR)=0.0 

35700 


8(NR)=-i.0+DX/0ELX 

35800 


C(iMR) = 1.0 

35900 


0(NR)=O.U 

36000 


00370 I = fiKP,NM 

36100 


A(r)=-ALA.ViS 

3fe200 


BU) = 1.0+2,{)*ALA/';S 

3b3uO 


C ( 1 )=-ALAmS 

3d4u0 

3/0 

i>tl)=T(I) 

36500 


A(^-i) = l .0 

36600 


BCi''i)=-l ,0 

36700 


CCrO=0,0 

36800 


DCiSiJsU.O 

36900 

380 

coin'icjUb; 

37000 


Du390 1=1, M 

37100 

390 

Oii(i)=U( fHCTH-TM) + IM 

37200 

V-. 

PHr,Ui'4uU,lDA(i),I = i,IM) 

373U0 

400 

F;JRviAtC5X,10(FlU.l,3X) ,//,5A,l0tFl0. 1 ,3X) ,//) 

374U0 


r i‘c;K2 = IlEH2i-l ,*iMiT=«lT+l 

37500 


call rHlDAG(A,o,C,D,N,XX3 

37600 


TlMi:I='r IME + OT 

37700 


0J420 i=l,N 

37800 


r(u=xxtij 



37QU0 

38000 

38100 

38200 

38300 

38400 

38500 

38600 

38700 

38800 

38900 

390u0 

39100 

39200 

39300 

39400 

39500 

396O0 

39700 

39800 

39900 

40000 

40100 

40200 

40300 

40400 

40500 

40600 

40700 

40800 

40900 

41000 

41100 

41700 

41800 

41900 

42000 

42100 


42200 


TA C 1) =T 1 1 ) * t TB-Tiw ) +TM 
PKINT430, CTACi 1,1=1, NR) 

FUHKATC5X,l0CFl0.i,3X),/) 

PRlAiT440,TCl ),T(NH5 ,TaO 
FJHmAT(5X, 'TEMPS= ' ,3CF10.2,3X) ,//) 

7.^ = (CONSORT CiMBP) )/(DXN-UF.LXN)+T(NHMJ/(DXNtOELXN) 

VfsCCOr.S'i'C lA(iJHP)-TM)/tUX-DEJuX)-C0NL*('fM-TA(wRM3)/tDX+UELX))/ 

9(DFNSi'ALMJ 

VA=(/C + VN 

0U = -CONi.*(TM-TA(NHM )/(uX+DFl,X) 
a3=:-C0N6«^(7’A(NHP)-TM)/(QX-De;LX) 

Ortt;L'r=0ENS*Ali.Vi*^A 
\/ = T(NHM)/(DXN + DtLXN) 

va=v^vcm 

DELX=VA#LVr 

Db'6X^'^ = DFL-A/AL. 

rS=i)Xr^'CTCl)-T(2))/DX+T(2) 

Oi6i=ul5TfDELX 

i'3A= C i , O^TA 1 1 ) -TA 1 2 ) ) /2 . 0 

IFtNI'r.UT.lNT2) GU lO 316 

Nrr=o 

PRIN ri65, flME,TSA,TA(l<) ,VA,OIST 
CONTINUE 

PRINTi65,TlME,TSA,TA(N3,VA,OIST,gL,QS,QMELT 

IFCT5.GE.1.) GO TO 530 

CONTINUE 

PRINT560 

F0RMA'T(5X,' VAporisation not started *,//) 

GO TO 520 
PRINr490 

FORMATClux, nvhole length melted ',//) 

CONTINUE 

CONTINUE 

PH1NT165,T1ME,T5A,TACn),VA,UIST 

PRlNT600,iTEtU,ITEF2 

F0R«AT(5X,' ITEkl= ',15,' 1TER2= ',15,//) 

PRINT540 

fokmaksx,' boiling starts ') 

STOP; END 



42400 

42500 

42600 

427UO 

42800 

429U0 

43000 

431O0 

43200 

43300 

43400 

43500 

43600 

43700 

438Q0 

43900 

44000 

44t00 

44200 

44300 

44400 
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Saoroutine to solve trlciiagondl matrix resulting from 
finite aitf, formulation, 

SuBHQUtlNfc: TR10AG(A,B,C,D,N,XX 
A(lu0o),aCl6003 ,ctio 
93AMAC1000) 

BfiTAd 

3AMACU=a(l)/B£TAt 1 J 
Dtj5 I™2rN 
J4=l-1 

HerACl3=BtI)-A(I)*C(I«)/B6:TA(iM) 

GAMA(I) = C0CiJ-Aa)4GAmAtI<4)l/»ETACIJ 
XXtR)=GAMA(N3 

N iJSN-l 
l)J10 1 = 1, Nf4 
K = H-1 

KP=K+1 , ^ 

XX(i<)=GAMAlK)-C(K)l'XXlKP)/BETA(K) 

RETURN 

E4D 


103 ,0(1000 ),XX(tOOOJ,OETA (1000 3, 



